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Abstract. Using extensive Molecular Dynamics simulations we study the 
behavior of polyclcctrolytes with hydrophobic side chains, which are known to 
form cylindrical micelles in aqueous solution. We investigate the stability of 
such bundles with respect to hydrophobicity, the strength of the electrostatic 
interaction, and the bundle size. We show that for the parameter range relevant 
for sulfonated poly-para-phenylenes (PPP) one finds a stable finite bundle size. 
In a more generic model we also show the influence of the length of the precursor 
oligomer on the stability of the bundles. We also point out that our model 
has close similarities to DNA solutions with added condensing agents, hinting 
to the possibility that the size of DNA aggregates is under certain circumstances 
thermodynamically limited. 
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1. Introduction 

Aggregates of charged semiflexible polymers have been experimentally observed in 
a variety of systems, where the most well known examples are probably toroidally 
shaped bundles in DNA solutions with added condensing agents, such as multivalent 
counterions PJGIE]- Another polyelectrolyte system of synthetic nature, consisting 
of poly(para-phenylenc) (PPP) oligomers @]03E], which are short rodlike charged 
objects, forms cylindrical micelles due to its hydrophobic side chains. This model 
can be well controlled chemically, and can be used as a synthetic model system to 
study nano-aggregation of charged filaments. An additional advantage of PPPs is that 
they allow one to tune experimentally the electrostatic and hydrophobic interactions 
separately. Although the physical origin of both aggregation mechanisms is different, 
in a first approximation they can be regarded as being the result of a short range 
attraction. This interaction is, in the case of DNA, the result of short range ionic 
correlations of the multivalent counterions 0, whereas in the case of the PPPs it is 
due to short range interactions of the hydrophobic side chains. The aggregate size in 
both cases is determined by the competition between the surface tension (due to short 
range attractions), the repulsive self-energy of the backbone charges, and the entropic 
degrees of freedom of the counterions. For DNA it has been argued |SJ E3 IU] 
that the observed finite size of the DNA aggregates is due to kinetic problems, 
i.e. electrostatic barrier formation, and not set by thermodynamic properties of the 
system. Experiments performed with viral DNA support this conjecture at least for 
the investigated systems ^2]- However, recently Henle and Pincus ^Sj have argued, 
using a charged rod model interacting via a short range attraction, that, depending 
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on the actual parameters of the system, either finite or infinite bundle sizes should 
be possible. Treating the systems as consisting of sticky charged rods brings up the 
analogy to the Rayleigh split of a charged hydrophobic droplet. This problem has 
been solved in a mean-field model already by Deserno He observes that the 

droplet size is always finite, if the counterions cannot penetrate into the droplet, but 
that allowing counterion penetration leads either to finite or infinite droplet sizes, 
depending on the parameters. There have been relatively few simulations on bundle 
formations, notable exception being the simulations by Stevens [151 llfij . showing the 
possibility that multivalent ions alone can lead to bundle formations, and the work 
of Borukhov et al. |17|. which investigates two-rod systems bundled via short range 
mobile linker interactions. In the present work we investigate the bundle formation 
for a system of charged semi-flexible polymers with short range interactions for two 
models, and demonstrate clearly the existence of finite size aggregates in both cases. 

2. Simulation method 

2.1. PPP 

The PPP oligomer is described with a bead spring model, whose mapping is shown 
in figure. ^ Each phenyl ring along the semi-flexible backbone is represented by a 
spherical bead. Two of them carry one negative unit charge e Q . The chain length 




Figure 1. Mapping of a PPP monomer onto the bead spring model. 



is N m — 61 beads representing 20 PPP monomers. The hydrophobic side chain 
is modeled as a flexible chain containing 3 beads. The counterions are simulated 
explicitly as charged spheres carrying one unit charge. The solvent is taken into 
account implicitly represented as a background with dielectric constant e. The 
simulations are performed in an NVT ensemble. The temperature is maintained via 
a Langevin thermostat. 

All particles interact via a repulsive Lennard-Jones potential with r cu t = 2 1 / 6 cr 
and clj = IksT (eq. and charged particles interact via an unscreened Coulomb 
potential (eq.[5J|. Bonds are realized with a FENE ^Hj potential with kp = 7 k bT and 
a cutoff Rf = 2cr (cq. 01). The PPP backbone is semi-flexible which is modeled with 
a bond angle potential with kg = 30fcsT (eq. 0}. The influence of the solvent on the 
hydrophobic side chains is taken care of with an attractive Lennard-Jones potential 
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between these monomers (eq. 0. Here the cutoff is r cut = 2.5a. 

C/Lj(r)=4e LJ ^) 12 -(^)^ for r < r cut (1) 

U C {r ij )=i B T^- (2) 

U F {r) = h F R 2 F In M - ^SSj ^ for r < R F (3) 

U e (r) = k g (l-cos0) (4) 



The parameters clj for the attractive Lennard- Jones interaction and the Bjerrum 
length £b = e\j{^KtkBT) have been varied. Due to the large kinetic barriers 
involved in the formation of a polyelectrolyte bundle, we study only its possible 
break up by starting the simulations with a preformed bundle in the middle of a 
spherical simulation cell. The cell radius is given by the particle density which is 
p = 1.0 x 10~ 4 <7~ 3 . The equilibration is done in a two step process. First we constrain 
the backbone monomers of the PPPs in space and let the hydrophobic hairs and 
counterions equilibrate for 100 000 MD steps. Then we release the backbone monomers 
and let the total system equilibrate. The bundle is simulated for 2 000 000 MD steps. 
If the bundle is stable over the whole simulation time it is called stable in this paper. 
We are well aware that this procedure only proves that the observed bundles are at 
least in a metastable state and cannot prove that we are really in the global minimum 
of the free energy. However, we only want to demonstrate the existence of finite 
aggregate sizes, and therefore our method suffices, see e.g. Ref. ^7] where the same 
simulation technique has been applied to the case of a two bundle system. 

2.2. Generic bundles 

Further simplification of this model can be achieved by taking the zero size limit of 
the hydrophobic side chains and making instead all backbone monomers hydrophobic, 
which is exactly the model Henle and Pincus suggested |13|. Each bead carries a 
positive unit charge. The beads interact with all other charges via an unscreened 
Coulomb interaction (eq. EJ). The beads are connected by a FENE potential (eq. 0J) 
and furthermore the oligomer is semi-flexible, where the stiffness of the chain can 
be tuned by a bond angle potential (eq. 0J. The hydrophobic interactions due to 
side chains are also represented via short range interactions of these beads, where the 
Lennard- Jones+Cosine potential (eq.EJ is used [d]. This potential enables the use of 
large e^j values, while the function and it's first derivative are continuous at the cut 
off radius. The cut off radius is chosen to be extremely small , r cut = 1.5a, so that 
the counterions can not penetrate into the bundle without breaking the short range 
interaction among the polyelectrolyte beads. All other interaction parameters are as 
described in section |2~T1 

Ucos{r) = (cos (ar 2 + 0) - l) for r min < r < r cut (5) 
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3. Results and discussions 

3.1. PPP 

For the PPP model we are interested in the influence of the hydrophobicity and the 
strength of the electrostatic interaction on the stability of the PPP micelles. We 
simulated PPP bundles of different aggregate sizes and varied systematically the 
parameters for the hydrophobicity e^j and the strength of the electrostatic interaction 
£ B - The simulated bundle sizes N p are 2, 3, 5, 8 and 10. In Figure |21 we show the 
stable and unstable regions in a phase diagram. 
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Figure 2. elj-^s phase diagram for bundles with N p =2, 3, 5 ,8 ,10. 

The phase boundary for bundles with different sizes are shown as lines connecting 
simulations where the bundles are marginally stable. Above the line bundles with that 
specific size are stable (actually, only at least meta-stable) , below unstable. 

Looking at the phase boundary for a given bundle size, e.g. N p — 5, with 
increasing £g the stability of the bundle first decreases. This means one needs a 
larger hydrophobicity to get a stable bundle. This can be explained by the increasing 
electrostatic repulsion of the charged backbones. At £b/ct ~ 2 the stability curve 
shows a maximum. Further increase of £3 leads to more stable bundles. This non- 
monotonic form resembles the behavior of the extension of flexible polyelectrolyte 
chains both in good |3U] and poor solvent |53J |^ . The increasing stability 

at high £b values is due to the increased density (condensation) of counterions close 
to the bundle which renormalize the effective charge and can also mediate attractive 
interactions at sufficiently large values of Is- 

In Figure |21 we can also compare the stability of bundles with different size. For 
£b = 0.5cr where the electrostatic interaction does not play a significant role, we see 
that for the investigated regime the range of e- values yielding stable bundles increases 
with increasing bundle size. This behavior changes drastically when one looks at £b 
values larger than 1.0a. Here we find at e-valucs between 1.7 and 1.9 that bundles 
which contain between 5 and 8 PPPs are more stable than larger or smaller N p values. 
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This implies also that one can find situations with a finite stable bundle size. Note 
that we observe this behavior in the relevant experimental region, namely Ib — 1.7 'o 
for PPP. We can explain this finding by an electrostatic instability in analogy to the 
Rayleigh instability of charged droplets [251 114). which has been shown to be valid 
also for linear charged polymers I^HJIJZ]. When one further increases Ib, the marginal 
stable bundle size decreases more and more, and the stability lines for all bundles 
containing more than 3 PPPs merge. We expect that for even larger £b values one 
will find an infinite bundle size. 

Finally we investigated the counterion distribution inside and in the close vicinity 
of a bundle. In Figure|3we show a snapshot of a polyelectrolyte bundle together with 
its cross section. Note that the counterions are able to penetrate also inside the bundle. 
This is important for the theoretical approach to the Rayleigh instability of charged 
droplets as it is introduced by Deserno [T4] , 




Figure 3. The snapshot of a polyelectrolyte bundle (left) together with a cross 
section (right) is shown for N p = 8, e^j = l.7kgT and £g = l.Tcr 



3.2. Generic bundles 

Also for the generic bundle model as discussed above we find the formation of finite 
size aggregates. In order to study the contribution of various effects (e.g. strength 
of Coulomb and Lennard-Jones interactions, stiffness of the chains) the stability of 
a bundle of certain size can be tested as follows. A cylindrical bundle is preformed 
by aligning the backbone of the polymers along the axis of the cylinder. The cross- 
section of the bundle is chosen as a close-packed structure on a triangular lattice 
to minimize the surface area. Unlike the PPP polymers in the previous section, 
where the hydrophobic side chains where explicit, in this more generic model the 
hydrophobic interaction is also associated with the backbone beads. Therefore the 
polymers are not packed like a classical micelle. The lattice spacing for this initial 
configuration is chosen as the minimum of the Lennard-Jones + Cosine potential. 
The counterions are placed in a cylinder enclosing the bundle, where the radius of 
this cylinder is large enough to provide sufficient freedom during the relaxation of the 
counterions. The equilibration is started with the counterions, which are set free to 
move within this cylinder surrounding the bundle. However, it is important to note 
that the counterions cannot penetrate into the bundle at this stage due to the close 
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packing of the polymers. After the counterions are equilibrated, the polymers are set 
free, and a brief equilibration run is performed where the polymers are also confined 
within a cylinder tightly enclosing the polymers, so that the bundle does not fall apart 
at this stage. Next the cylinder constraints are removed, and the system is integrated 
for 2.0 x 10 6 time steps. During the simulation, depending on the relative value of 
the interaction parameters, the bundle either remains as a single aggregate or splits 
up into smaller aggregates. 

In Figure0the biggest stable bundle size obtained from the simulation is shown as 
a function of the number of polymers in the initial bundle. The polymers are iV m =10 
beads long, and e^j = A.OksT and Ib = 2. Oct. The stiffness of the chains is kept 
fixed (kg = lOfcsT). Up to six polymers per startup bundle, the structure does not 
disintegrate and the bundle remains as a single entity. However, for initial aggregates 
with more than N p =6 polymers, the aggregate splits into two or more bundles. In 
Figure the biggest size among these bundles is plotted. The stable equilibrium 
bundle size for this set of parameters is N p « 5. However, due to the finite system size 
the equilibrium aggregate size obtained via this method must be taken cautiously. 
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Figure 4. The biggest stable aggregate size (JV™ a:c ) as a function of number of 
polymers (N p ) in the startup bundle. 

Another approach to obtain the average bundle size is to look at the internal 
energy of a bundle with counterions. The bundle is formed as in the previous case; 
however, in this simulation the bundle is constraint to remain as a single aggregate. 
In other words, the polymers are fixed throughout the simulation, and only the 
counterions are free to move. The average internal energy of the system as a function 
of the number of polymers per aggregate is plotted in Figure[S]for the same interaction 
parameters as in the previous case. In these constrained bundle simulations a minimum 
in the internal energy of the system is observed for N p =3. Initially as the number 
of polymers per bundle increases, the internal energy of the system decreases as a 
result of decreasing surface energy. However, above the average aggregate size the 
electrostatic energy of the system increases dramatically. The energy distribution for 
bundles smaller than N p =7 polymers is relatively flat. The sudden increase in the 
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internal energy for the N p =7 polymer aggregate is due to the complete isolation of 
the central polymer from the surrounding counterions. This is an artifact of fixing 
the polymers during the simulation, ft is important to note that only the internal 
energy of the system is taken into account in this analysis. However, the entropy of 
the counterions could also play a dominant role in determining the average aggregate 
size for these polyelectrolytes with hydrophobic interactions. The slight difference 
in the equilibrium aggregate size obtained with these two methods presented above 
can be explained by the entropic contribution which is neglected in the fixed bundle 
simulations. 
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Figure 5. The average internal energy per polymer as a function of number of 
polymers (N p ) in the startup bundle. 

The length and charge density of the precursor polymers effect the fraction of 
condensed ions. Even though Poisson-Boltzmann theory yields a good approximation 
for the case of infinitely long polymers for finite length rods there are no straight- 
forward solutions. In order to gain further inside into the finite length effects, we have 
performed bundle simulations for lengths ranging from iV m =10 to iV m =60 monomers. 
Snapshots from the equilibrated systems of iV m =20, N m =40, and N m —60 are shown 
in FigureEl For bundles with less than N p =6 polymers all simulated polymer lengths 
yield stable aggregates, no splitting is observed. On the other hand, for bundles of 6 
polymers even though the polymers of length N m =10 and iV m =20 monomer length 
provide stable aggregates, for longer polymer lengths the bundle splits up. 

The internal energy and the contributions from Coulombic and Lennard- Jones 
potentials per particle for bundles of N p =6 polymers with polymer length of N m =60 
monomers is shown in Figure \7\ Upon splitting, the internal energy of the bundle 
increases by almost one UbT per particle. The split takes place rather abruptly in 
a late stage of the simulation, which suggests the presence of a high energy barrier 
for the split. The split of the bundle into two small bundles is unfavorable in terms 
of the short range Lennard- Jones interactions since the surface energy increases. On 
the other hand for the Coulombic interactions the answer is more complicated. The 
electrostatic self energy of the bundle decreases due to split, since the like-charged 
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Figure 6. Snapshots from simulations of a 6 polymer bundle with polymer 
lengths N m =20, A r m =40, and ./V m =60 from left to right, respectively. 



polymers repel each other. However, since the split-bundles attract less covmtcrions 
the favorable electrostatic interactions among the polymers and counterions is also 
lost. For N p =G polymer bundle upon splitting the Coulombic contribution to the 
internal energy decreases. 




r, time 



Figure 7. The total internal energy (middle), Coulomb (top) and Lennard- 
Jones (bottom) contributions per monomer for bundles of 6 polymers during the 
simulation. 

The increase in the internal energy of the system can be advantageous for the 
system only if the release of the counterions sufficiently increases the entropy to 
compensate for the increase in internal energy. In Figure [S] the integrated counterion 
distribution for the N p —6 polymer bundle is given for five time intervals during the 
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simulation. For each counterion the distance to the closest monomer is chosen for 
the radial distribution function. Prior to the split (s0-s3) the counterion distribution 
does not change. Upon splitting up into two bundles the fraction of ions close to the 
bundles decreases dramatically. Therefore we can conclude that the splitting of these 
long polymer bundles are driven by the entropy of the released counterions. 




1 10 100 

r, distance from bundle 



Figure 8. Integrated counterion distribution as a function of radial distance 
from the 6 polymer bundle. The distributions are calculated for five different 
time periods during the simulation. Time periods from top curve to the bottom 
one: sO (0-399) - s4 (1600-1999). Time period s4 corresponds to the distribution 
after the bundle splits up. 



4. Conclusions 

We have shown that charged semi-flexible polyelectrolytes interacting via two kinds 
of hydrophobic short range interactions have parameter regions where finite aggregate 
sizes exist. For large values of Is the bundle size increases, and in the limit of 
large Is we observe the expected trend to form infinite bundles due to counterion 
crystallization. Furthermore, we could demonstrate that the release of counterions 
during a bundle split can considerably lower the free energy of the total system. 

The stability curve for different bundle sizes as a function of I b shows similarities 
to the Rayleigh instability, i.e., for an increase in Is one needs an increased clj. 
Another observation is that the stability line varies non-monotonically with Bjerrum 
length, which resembles the non-monotonic extension behavior linear polyelectrolytes 
show if the Bjerrum length is increased. We relate this to the fact that the electrostatic 
self-energy of the bundle does not increase with lg after 1b ~ 2er due to counterion- 
backbone ion correlations. The degree of polymerization shows also an effect on bundle 
size, which has been observed in experiments @J |S] on PPPs as well; however, our 
limited data does not allow us to draw definite conclusions about the underlying 
mechanism. 

One important point to note is the similarities of our model with the observed 
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trends in DNA condensation experiments which also find finite aggregate sizes. The 
applicability of our model to those experiments rests basically on two assumptions. 
First, we assume that the strength of the short range attraction does not change 
with external parameters, and second, presently we have not gauged the interaction 
strength to that originating from multivalent counterion interactions. This important 
piece of work is left for future investigations. However, it appears plausible that at 
least for some parameter regions of biological charged polymers both assumptions are 
justified, so that our results can be used to provide a mechanism for finite bundle 
sizes. 

Our results also demonstrate that the mean-field theory of Deserno ^31 f° r 
a charged hydrophobic droplet looks qualitatively correct. Extending the Deserno 
analysis to charged sticky rods does not change the qualitative behavior |29| . 
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